rm(list=ls())

##############################################################################
# FIGURE 1 SCRIPT
##############################################################################

# load packages
#library(ggplot2)
#library(glm.predict)

# load data 
load("data/candidate_GCI.RData")

##############
# RUN MODELS #
##############

LM1 <- lm(votes ~ ranking + factor(region), data = candidate_GCI)
#summary(LM1)
LM2 <- lm(votes ~ ranking + factor(region) + factor(elect_year), data = candidate_GCI)
#summary(LM2)
LM3 <- lm(votes ~ ranking + factor(region) + log_gdp_percap + 
            log_pop + idealpoint_US + idealpoint_RU + nato_bin + 
            eu_bin + g77_bin + oic_bin, data = candidate_GCI)
#summary(LM3)

###############
# CREATE PLOT #
###############

# create case matrix
intercept <- rep(1, 140)
rank <- c(17:156)
regB <- rep(0, 140)
regC <- rep(0, 140)
regD <- rep(1, 140)
regE <- rep(0, 140)

cases <- cbind(intercept, rank, regB, regC, regD, regE)

# check range of rakings in RegionD
regionD <- subset(candidate_GCI, region == "D")
# 17 to 156

# create pred values for LM1 
pred_vals_matrix <- matrix(nrow = 140, ncol = 3)

for(i in 1:140){
  predvals <- glm.predict::basepredict(LM1, 
                                       cases[i,], 
                                       sim.count=1000, 
                                       conf.int=0.95, 
                                       set.seed=123,
                                       type = c("simulation"), 
                                       summary = TRUE)
  pred_vals_matrix[i, ] <- predvals
}

# create pred val plot for LM1 
plot_data3 <- as.data.frame(cbind(pred_vals_matrix, cases))
names(plot_data3)[1:3] <- c("estimate", "lower", "upper")

gg_CII1 <- ggplot(plot_data3, aes(x = rank, y = estimate)) + 
  geom_line() + 
  geom_ribbon(aes(ymin = lower, ymax = upper), 
              alpha = 0.5) + 
  labs(
    y = "Predicted Number of Votes", x = "GCI Ranking", 
    caption = "Typical case: Region = Africa") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5), legend.position = 'none') + 
  theme_bw() + geom_rug(data = subset(candidate_GCI, region == "D"), aes(x = ranking),inherit.aes = FALSE, color = "black") 

gg_CII1

# save pdf
ggsave(gg_CII1, filename = "output/fig1.pdf")

# script complete message 
print("fig1 complete")

